      SUBROUTINE FCDSK(RHO,TH,PH,EPS,FV)
      PARAMETER MAX=25
      DIMENSION FO(3)
      DOUBLE PRECISION FV(3),XMU,RC
      DOUBLE PRECISION RHO,TH,PH,EPS,FO,CO,THRAD,PHRAD,RCOS,RP(3),PLCO,T
      DOUBLE PRECISION ALI(3,3),PHII(3,3),AI(3,3),PLI(3,3),RI(3),THETA,L
      DOUBLE PRECISION DILAT,DLCO
      DOUBLE PRECISION DELGAM,GAMMA1,RM2,FATSO,LIBBER,DELGM1,       N,M
      DOUBLE PRECISION JONE,JTWO,ION1,ITW2,ITH3EE,GRHO,GTHETA,G,FORCE(3)
      DOUBLE PRECISION KONE,     IONE,ITWO,ITHREE
      DOUBLE PRECISION ZI1(MAX),ZI2(MAX),ZI3(MAX) ,FACTOR
      DOUBLE PRECISION E,K2,EONE
      DOUBLE PRECISION COM(1)
      DOUBLE PRECISION T2P1
      DATA CO /.1745392925199430D-1/
      COMMON RC,XMU,COM
      COMMON /CRVDK/AMASS,IDSWTH,NUMRN(1)
      COMMON /PRAM/ NMASS,NHARM,NDISK,NCDSK
      LOGICAL IDSWTH
      IDING=NMASS*4+2*NHARM*NHARM+NHARM+5*NDISK
      INTERG=12
      G=XMU
      E=0.0
      CRIT=1.D-2
      THRAD=TH*CO
      PHRAD=PH*CO
      RCOS=RHO*DCOS(THRAD)
      RP(1)=RCOS*DCOS(PHRAD)
      RP(2)=RCOS*DSIN(PHRAD)
      RP(3)=RHO*DSIN(THRAD)
      DO 10 I=1,3
      FV(I)=0.D0
      DO 20 J=1,3
      ALI(I,J)=0.D0
      PHII(I,J)=0.D0
      AI(I,J)=0.D0
   20 CONTINUE
   10 CONTINUE
      PHII(2,2)=1.D0
      ALI(3,3)=1.D0
      IRNUM=0
      IDING=IDING-3
      DO 100 I=1,NCDSK
      IDING=IDING+4+IRNUM
      ILAT=IDING
      ILON=IDING+1
      IRADC=IDING+2
      IRNGD=IDING+3
      DLCO=COM(ILON)*CO
      ALI(1,1)=DCOS(DLCO)
      ALI(2,2)=ALI(1,1)
      ALI(1,2)=DSIN(DLCO)
      ALI(2,1)=-ALI(1,2)
      PLCO=COM(ILAT)*CO
      PHII(1,1)=DSIN(PLCO)
      PHII(3,3)=PHII(1,1)
      PHII(3,1)=DCOS(PLCO)
      PHII(1,3)=-PHII(3,1)
      CALL MATMUL(PHII,ALI,3,3,3,PLI)
      CALL MATVEC(PLI,RP,3,3,RI)
      RM2=RI(1)*RI(1)+RI(2)*RI(2)+RI(3)*RI(3)
      RM2=DSQRT(RM2)
      T=COM(IRADC)/RM2
      T3=T*T*T
      T2P1=T*T+1
      THETA=3.1415926536D0 /2.0-DASIN(RI(3)/RM2)
      LIBBER=2.*T*DSIN(THETA)
      IONE=0.D0
      ITWO=0.0D0
      ITHREE=0.0D0
      FATSO=2.*T*DCOS(THETA)
      INTERI=2*INTERG
      INTERA=INTERI-1
      INTERB=INTERG-1
      INTERC=INTERI+1
      IRNUM=NUMRN(I)
      DELGAM=COM(IRNGD)*CO/(IRNUM*INTERI)
      GAMMA1=DELGAM*(IRNUM*INTERI+1)
      DO 200 JACKOF=1,IRNUM
      IDENSE=IRNGD+JACKOF
      DENSE=COM(IDENSE)
      IF (DENSE .EQ. 0.D0) GO TO 200
      IF(IDSWTH) DENSE=DENSE*1.0D10/AMASS
      DO 250 LDIVRN=1,INTERC
      GAMMA1=GAMMA1-DELGAM
      N=T2P1-FATSO*DCOS(GAMMA1)
      L=LIBBER*DSIN(GAMMA1)/N
      M=4./(N*DSQRT(N))
      K2=2.*L/(1.+L)
      CALL DCOKE(K2,KONE,EONE)
      JONE=M/(DSQRT(1.+L)*(1.-L))*EONE
      JTWO=0.D0
      IF(L.EQ.0.D0) GO TO 40
      JTWO=-(M*((1-L)*KONE-EONE))/(L*DSQRT(1.+L)*(1-L))*DSIN(GAMMA1)
      IF(K2 .LT. 2.0D-17) JTWO=0.D0
   40 ZI1(LDIVRN)=DENSE*DCOS(GAMMA1)*DSIN(GAMMA1)*JONE
      ZI2(LDIVRN)=DENSE*DSIN(GAMMA1)*JONE
      ZI3(LDIVRN)=DENSE*DSIN(GAMMA1)*JTWO
      GO TO 250
  250 CONTINUE
      DELGM1=DELGAM*INTERI
      ION1=0.D0
      ITW2=0.D0
      ITH3EE=0.D0
      DO 50 IWONG=1,INTERB
      ION1=ION1+4.*ZI1(2*IWONG)+2.*ZI1(2*IWONG+1)
      ITW2=ITW2+4.*ZI2(2*IWONG)+2.*ZI2(2*IWONG+1)
      ITH3EE=ITH3EE+4.*ZI3(2*IWONG)+2.*ZI3(2*IWONG+1)
   50 CONTINUE
      ION1=ION1+4.*ZI1(INTERI  )+ZI1(1)+ZI1(INTERC)
      ITW2=ITW2+4.*ZI2(INTERI  )+ZI2(1)+ZI2(INTERC)
      ITH3EE=ITH3EE+4.*ZI3(INTERI  )+ZI3(1)+ZI3(INTERC)
      FACTOR=DELGM1/(6.*INTERG)
      ION1=ION1*FACTOR
      ITW2=ITW2*FACTOR
      ITH3EE=ITH3EE*FACTOR
      IONE=IONE+ION1
      ITWO=ITWO+ITW2
      ITHREE=ITHREE+ITH3EE
      GAMMA1=GAMMA1+DELGAM
  200 CONTINUE
      GRHO=-G*T*T*T*(1./T*ITWO-DCOS(THETA)*IONE-DSIN(THETA)*ITHREE)
      GTHETA=-G*T*T*T*(DSIN(THETA)*IONE-DCOS(THETA)*ITHREE)
      DILAT=DSQRT(RI(1)*RI(1)+RI(2)*RI(2))
      FORCE(1)=1./DILAT*(GTHETA*DCOS(THETA)+GRHO*DSIN(THETA))
      FORCE(2)=RI(2)*FORCE(1)
      FORCE(1)=RI(1)*FORCE(1)
      FORCE(3)=-GTHETA*DSIN(THETA)+GRHO*DCOS(THETA)
      DO 300 IMESS=1,3
      DO 310 J=1,3
      AI(IMESS,J)=PLI(J,IMESS)
  310 CONTINUE
  300 CONTINUE
      CALL MATVEC(AI,FORCE,3,3,FO)
      DO 340 JSCUM=1,3
      FV(JSCUM)=FV(JSCUM)+FO(JSCUM)
  340 CONTINUE
  100 CONTINUE
      RETURN
      END
